#====# Appendix M: Heterogeneous effects #====#

# Load libraries and set defaults ----
library(bizdays)
library(fixest)
library(interflex)
library(marginaleffects)
library(modelsummary)
library(tinytable)
library(tidyverse)
library(tidylog, warn.conflicts = FALSE)
source("aux/plot_theme.R")
source("aux/event_function.R")

# prepare business calendar:
business_calendar <- create.calendar('biz_calendar', weekdays = c('saturday','sunday'))

# Import data ----
stocks <- read_rds("data_out/stocks_analysis.rds") # main analysis dataset

# Table M.1: Heterogeneous effect on past FCPA targets by whether they suffered enforcement action or only investigation ----
est <- return_daily_avg(data = stocks %>%
                          filter(FCPA_sample == 1), 
                        abn_ret = "abn_chg", cum_abn_ret = "car", stratum = "enforcement_action") %>%
  mutate(date = format(date, format = "%a, %b %d %Y")) %>%
  rename("term" = "date")

mod_ar1 <- list(tidy = est %>%
                  filter(enforcement_action == "1" &
                           dv == "ar") %>%
                  select(-enforcement_action, -dv),
                glance = est %>%
                  filter(enforcement_action == "1" &
                           dv == "ar") %>%
                  filter(term == "Mon, Feb 10 2025") %>%
                  select(N) %>%
                  rename("N of firms" = "N")) %>%
  `class<-`("modelsummary_list")

mod_car1 <- list(tidy = est %>%
                   filter(enforcement_action == "1" &
                            dv == "car") %>%
                   select(-enforcement_action, -dv),
                 glance = est %>%
                   filter(enforcement_action == "1" &
                            dv == "car") %>%
                   filter(term == "Mon, Feb 10 2025") %>%
                   select(N) %>%
                   rename("N of firms" = "N")) %>%
  `class<-`("modelsummary_list")

mod_ar0 <- list(tidy = est %>%
                  filter(enforcement_action == "0" &
                           dv == "ar") %>%
                  select(-enforcement_action, -dv),
                glance = est %>%
                  filter(enforcement_action == "0" &
                           dv == "ar") %>%
                  filter(term == "Mon, Feb 10 2025") %>%
                  select(N) %>%
                  rename("N of firms" = "N")) %>%
  `class<-`("modelsummary_list")

mod_car0 <- list(tidy = est %>%
                   filter(enforcement_action == "0" &
                            dv == "car") %>%
                   select(-enforcement_action, -dv),
                 glance = est %>%
                   filter(enforcement_action == "0" &
                            dv == "car") %>%
                   filter(term == "Mon, Feb 10 2025") %>%
                   select(N) %>%
                   rename("N of firms" = "N")) %>%
  `class<-`("modelsummary_list")

mod_ar_diff <- list(tidy = est %>%
                      filter(enforcement_action == "diff" &
                               dv == "ar") %>%
                      select(-enforcement_action, -dv),
                    glance = est %>%
                      filter(enforcement_action == "diff" &
                               dv == "ar") %>%
                      filter(term == "Mon, Feb 10 2025") %>%
                      select(N) %>%
                      rename("N of firms" = "N")) %>%
  `class<-`("modelsummary_list")

mod_car_diff <- list(tidy = est %>%
                       filter(enforcement_action == "diff" &
                                dv == "car") %>%
                       select(-enforcement_action, -dv),
                     glance = est %>%
                       filter(enforcement_action == "diff" &
                                dv == "car") %>%
                       filter(term == "Mon, Feb 10 2025") %>%
                       select(N) %>%
                       rename("N of firms" = "N")) %>%
  `class<-`("modelsummary_list")

modelsummary(list("(1) \\textsc{ar}" = mod_ar1,
                  "(2) \\textsc{car}" = mod_car1,
                  "(3) \\textsc{ar}" = mod_ar0,
                  "(4) \\textsc{car}" = mod_car0,
                  "(5) \\textsc{ar}" = mod_ar_diff,
                  "(6) \\textsc{car}" = mod_car_diff),
             # statistic = "[{conf.low}, {conf.high}]",
             notes = paste0("Average \\textsc{ar} and \\textsc{car} to past FCPA targets that have been under enforcement action or only investigated, ",
                            "per day. Standard errors of the mean reported in parentheses. P-values from a two-tailed test ",
                            "of difference from zero for the average against a standard normal distribution. ",
                            "Estimation window starts 30 days and ends 5 days before FCPA Executive Order. Market models estimated using ",
                            "the LASSO and individual S\\&P 500 constituents as predictors, selected using 15-fold cross validation. ",
                            "Columns 5 and 6 report the difference in means, respectively, between averages in columns 1 and 3, and those in columns 2 and 4.",
                            collapse = ""),
             title = "Heterogeneous effect on past FCPA targets by whether they suffered enforcement action or only investigation \\label{tab:het_enf_inv}",
             stars = c("*" = 0.05), 
             escape = FALSE) %>%
  group_tt(j = list("FCPA Enforcement" = 2:3,
                    "FCPA investigation" = 4:5,
                    "Difference-in-means" = 6:7)) %>%
  group_tt(i = list("Pre-event:" = 1,
                    "Post-event:" = 11)) %>%
  style_tt(i = 21, line_color = "white", line_width = 0.1, line = "t") %>%
  style_tt(i = 22, line_color = "black", line_width = 0.05, line = "b") %>%
  theme_tt("resize", width = .9) %>%
  theme_tt("placement", latex_float = "!htbp") %>%
  save_tt("tables/table_M1.html", overwrite = TRUE)

# Figure M.1: Average AR of past FCPA target firms by the decade in which the firm experienced its last FCPA action ----
summary(stocks$latest_year)

stocks %>%
  filter(latest_year < 2000) %>%
  filter(date == as.Date("2025-02-10")) %>%
  select(conm, ticker_symbol, latest_year, abn_chg)
# only 5 firms. We can group them with the 2000s

# what has missing value?
stocks %>%
  filter(FCPA_sample == 1) %>%
  filter(is.na(latest_year)) %>%
  filter(date == as.Date("2025-02-10")) %>%
  select(conm, ticker_symbol, latest_year, abn_chg)
# only CHEMRING, which is genuinely not reporting information on the FCPA Clearinghouse

dat <- stocks %>%
  filter(date == as.Date("2025-02-10")) %>%
  filter(FCPA_sample == 1) %>%
  mutate(decade = case_when(latest_year < 2010 ~ "action year < 2010",
                            latest_year < 2020 ~ "2010 ≤ action year < 2020",
                            latest_year <= 2025 ~ "action year ≥ 2020") %>%
           fct_relevel("action year < 2010", "2010 ≤ action year < 2020", "action year ≥ 2020"))

out <- map(.x = unique(dat$decade[!is.na(dat$decade)]),
           .f = function(x){
             d <- dat %>%
               filter(decade == x)
             
             t <- t.test(d$abn_chg)
             
             out <- data.frame(estimate = t$estimate,
                               std.err = t$stderr,
                               p.value = t$p.value,
                               N = sum(!is.na(d$abn_chg)),
                               decade = x) %>%
               `rownames<-`(NULL)
             return(out)}) %>%
  bind_rows() %>%
  mutate(conf.low = estimate - 1.96*std.err,
         conf.high = estimate + 1.96*std.err,
         conf.low.diff = estimate - qnorm(1-(1-0.834)/2)*std.err,
         conf.high.diff = estimate + qnorm(1-(1-0.834)/2)*std.err,
         lab = paste0("N = ", N)) %>%
  arrange(decade) 

# calculate p-values of pairwise difference-in-means tests:
p1 <- t.test(dat$abn_chg[dat$decade == "action year < 2010"],
             dat$abn_chg[dat$decade == "2010 ≤ action year < 2020"])$p.value

p2 <- t.test(dat$abn_chg[dat$decade == "2010 ≤ action year < 2020"],
             dat$abn_chg[dat$decade == "action year ≥ 2020"])$p.value

p3 <- t.test(dat$abn_chg[dat$decade == "action year < 2010"],
             dat$abn_chg[dat$decade == "action year ≥ 2020"])$p.value

out %>%
  mutate(decade = str_remove(decade, "action ") %>%
           fct_relevel("year < 2010", "2010 ≤ year < 2020")) %>%
  ggplot(aes(x = decade, y = estimate)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
  geom_linerange(aes(ymin = conf.low.diff, ymax = conf.high.diff), linewidth = 1) +
  geom_hline(yintercept = 0, linewidth = 0.1) +
  geom_text(aes(label = lab), hjust = -0.5) +
  # first horizontal line with vertical dents and p-value label:
  annotate(geom = "segment", 
           x = 1, xend = 2, y = 1.5, linewidth = 0.2) +
  geom_segment(data = data.frame(x = c(1, 2),
                                 xend = c(1, 2),
                                 y = c(1.5, 1.5),
                                 yend = c(1.48, 1.48)),
               mapping = aes(x = x, xend = xend, y = y, yend = yend), linewidth = 0.2) +
  annotate(geom = "text",
           x = 1.5, y = 1.56, label = paste0("p = ", sprintf("%.2f", p1))) +
  # second horizontal line with vertical dents and p-value label:
  annotate(geom = "segment",
           x = 2, xend = 3, y = 1.6, linewidth = 0.2) +
  geom_segment(data = data.frame(x = c(2, 3),
                                 xend = c(2, 3),
                                 y = c(1.6, 1.6),
                                 yend = c(1.58, 1.58)),
               mapping = aes(x = x, xend = xend, y = y, yend = yend), linewidth = 0.2) +
  annotate(geom = "text",
           x = 2.5, y = 1.66, label = paste0("p = ", sprintf("%.2f", p2))) +
  # third horizontal line with vertical dents and p-value label:
  annotate(geom = "segment",
           x = 1, xend = 3, y = -0.3, linewidth = 0.2) +
  geom_segment(data = data.frame(x = c(1, 3),
                                 xend = c(1, 3),
                                 y = c(-0.3, -0.3),
                                 yend = c(-0.28, -0.28)),
               mapping = aes(x = x, xend = xend, y = y, yend = yend), linewidth = 0.2) +
  annotate(geom = "text",
           x = 2, y = -0.35, label = paste0("p = ", sprintf("%.2f", p3))) +
  xlab("Year of last experienced FCPA action (enforcement or investigation)") + 
  ylab("Average abnormal returns on Executive Order day") +
  scale_y_continuous(labels = ~paste0(.x, "%"))
ggsave("plots/figure_M1.pdf", width = 7, height = 5)

# Figure M.2: verage AR of past FCPA target firms by five-year period in which the firm experienced its last FCPA action ----
dat <- stocks %>%
  filter(date == as.Date("2025-02-10")) %>%
  filter(FCPA_sample == 1) %>%
  filter(latest_year >= 2000) %>%
  mutate(period = case_when(latest_year < 2005 ~ "action year < 2005",
                            latest_year < 2010 ~ "2005 ≤ action year < 2010",
                            latest_year < 2015 ~ "2010 ≤ action year < 2015",
                            latest_year < 2020 ~ "2015 ≤ action year < 2020",
                            latest_year <= 2025 ~ "action year ≥ 2020") %>%
           fct_relevel("action year < 2005", "2005 ≤ action year < 2010",
                       "2010 ≤ action year < 2015", "2015 ≤ action year < 2020",
                       "action year ≥ 2020"))

out <- map(.x = unique(dat$period[!is.na(dat$period)]),
           .f = function(x){
             d <- dat %>%
               filter(period == x)
             
             t <- t.test(d$abn_chg)
             
             out <- data.frame(estimate = t$estimate,
                               std.err = t$stderr,
                               p.value = t$p.value,
                               N = sum(!is.na(d$abn_chg)),
                               period = x) %>%
               `rownames<-`(NULL)
             return(out)}) %>%
  bind_rows() %>%
  mutate(conf.low = estimate - 1.96*std.err,
         conf.high = estimate + 1.96*std.err,
         conf.low.diff = estimate - qnorm(1-(1-0.834)/2)*std.err,
         conf.high.diff = estimate + qnorm(1-(1-0.834)/2)*std.err,
         lab = paste0("N = ", N)) %>%
  arrange(period) 

out %>%
  mutate(period = str_remove(period, "action ") %>%
           fct_relevel("year < 2005", "2005 ≤ year < 2010",
                       "2010 ≤ year < 2015", "2015 ≤ year < 2020")) %>%
  ggplot(aes(x = period, y = estimate)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
  geom_linerange(aes(ymin = conf.low.diff, ymax = conf.high.diff), linewidth = 1) +
  geom_hline(yintercept = 0, linewidth = 0.1) +
  geom_text(aes(label = lab), hjust = -0.5) +
  xlab("Year of last experienced FCPA action (enforcement or investigation)") + 
  ylab("Average abnormal returns on Executive Order day") +
  scale_y_continuous(labels = ~paste0(.x, "%"))
ggsave("plots/figure_M2.pdf", width = 8, height = 5)

# Table M.2: Heterogeneous effect of Executive Order by year of last FCPA action----
dat <- stocks %>%
  filter(date >= offset(event, -5, business_calendar) &
           date <= offset(event, +5, business_calendar)) %>%
  mutate(event_day = as.numeric(date == event)) %>%
  filter(FCPA_sample == 1) 

mod <- feols(abn_chg ~ latest_year*event_day,
             data = dat,
             cluster = ~gvkey)

rows <- tibble::tribble(
  ~term, ~mod,
  "Std.Errors", "by: firm")
attr(rows, "position") <- c(9)

modelsummary(list("(1)" = mod),
             stars = c("*" = 0.05),
             escape = FALSE,
             coef_map = c("latest_year:event_day" = "Executive Order $\\times$ Year of FCPA action",
                          "event_day" = "Executive Order",
                          "latest_year" = "Year of FCPA action",
                          "(Intercept)" = "(Intercept)"),
             gof_omit = "R2 Within|AIC|BIC|RMSE|FE|Std.Errors",
             notes = paste0("Linear models of \\textsc{ar} to past FCPA targets per day. ",
                            "Standard errors clustered by firm in parentheses. ",
                            "Estimation window starts 30 days and ends 5 days before FCPA Executive Order. Market models estimated using ",
                            "the LASSO and individual S\\&P 500 constituents as predictors, selected using 15-fold cross validation. ",
                            collapse = ""),
             title = "Heterogeneous effect of Executive Order by year of last FCPA action \\label{tab:het_year}",
             add_rows = rows) %>%
  group_tt(j = list("Past FCPA targets" = 2)) %>%
  theme_tt("resize", width = .9) %>%
  theme_tt("placement", latex_float = "!htbp") %>%
  save_tt("tables/table_M2.html", overwrite = TRUE)

# Figure M.3: Marginal effect of the FCPA executive order on past FCPA targets’ AR conditional on the year of last FCPA action ----
inter <- interflex(estimator = "binning",
                   data = as.data.frame(dat),
                   Y = "abn_chg", D = "event_day", X = "latest_year",
                   na.rm = TRUE, vcov.type = "cluster", cl = "gvkey",
                   figure = FALSE)

shift <- 1
avg_comparisons(mod, variables = "event_day", by = "latest_year", 
                type = "response",
                newdata = datagrid(latest_year = unique(stocks$latest_year[!is.na(stocks$latest_year)]),
                                   event_day = c(0, 1))) %>%
  ggplot() +
  geom_line(aes(x = latest_year, y = estimate+shift)) +
  geom_ribbon(aes(x = latest_year, ymin = conf.low+shift, ymax = conf.high+shift),
              alpha = .3) +
  geom_pointrange(data = inter$est.bin$`1`,
                  mapping = aes(x = x0, y = coef+shift,
                                ymin = CI.lower+shift, ymax = CI.upper+shift)) +
  geom_hline(yintercept = 0+shift, linewidth = 0.1) +
  geom_histogram(data = dat,
                 mapping = aes(x = latest_year, y = 3*after_stat(count / sum(count)))) +
  scale_y_continuous("Predicted effect of Executive Order on abnormal returns",
                     labels = ~paste0(.x-shift, "%")) +
  xlab("Year of last experienced FCPA action (enforcement or investigation)")
ggsave("plots/figure_M3.pdf", width = 6, height = 5)

# Table M.3: Heterogeneous effect on past FCPA targets by whether they changed name since FCPA action ----
est <- return_daily_avg(data = stocks %>%
                          filter(FCPA_sample == 1), 
                        abn_ret = "abn_chg", cum_abn_ret = "car", stratum = "ticker_change") %>%
  mutate(date = format(date, format = "%a, %b %d %Y")) %>%
  rename("term" = "date")

mod_ar1 <- list(tidy = est %>%
                  filter(ticker_change == "1" &
                           dv == "ar") %>%
                  select(-ticker_change, -dv),
                glance = est %>%
                  filter(ticker_change == "1" &
                           dv == "ar") %>%
                  filter(term == "Mon, Feb 10 2025") %>%
                  select(N) %>%
                  rename("N of firms" = "N")) %>%
  `class<-`("modelsummary_list")

mod_car1 <- list(tidy = est %>%
                   filter(ticker_change == "1" &
                            dv == "car") %>%
                   select(-ticker_change, -dv),
                 glance = est %>%
                   filter(ticker_change == "1" &
                            dv == "car") %>%
                   filter(term == "Mon, Feb 10 2025") %>%
                   select(N) %>%
                   rename("N of firms" = "N")) %>%
  `class<-`("modelsummary_list")

mod_ar0 <- list(tidy = est %>%
                  filter(ticker_change == "0" &
                           dv == "ar") %>%
                  select(-ticker_change, -dv),
                glance = est %>%
                  filter(ticker_change == "0" &
                           dv == "ar") %>%
                  filter(term == "Mon, Feb 10 2025") %>%
                  select(N) %>%
                  rename("N of firms" = "N")) %>%
  `class<-`("modelsummary_list")

mod_car0 <- list(tidy = est %>%
                   filter(ticker_change == "0" &
                            dv == "car") %>%
                   select(-ticker_change, -dv),
                 glance = est %>%
                   filter(ticker_change == "0" &
                            dv == "car") %>%
                   filter(term == "Mon, Feb 10 2025") %>%
                   select(N) %>%
                   rename("N of firms" = "N")) %>%
  `class<-`("modelsummary_list")

mod_ar_diff <- list(tidy = est %>%
                      filter(ticker_change == "diff" &
                               dv == "ar") %>%
                      select(-ticker_change, -dv),
                    glance = est %>%
                      filter(ticker_change == "diff" &
                               dv == "ar") %>%
                      filter(term == "Mon, Feb 10 2025") %>%
                      select(N) %>%
                      rename("N of firms" = "N")) %>%
  `class<-`("modelsummary_list")

mod_car_diff <- list(tidy = est %>%
                       filter(ticker_change == "diff" &
                                dv == "car") %>%
                       select(-ticker_change, -dv),
                     glance = est %>%
                       filter(ticker_change == "diff" &
                                dv == "car") %>%
                       filter(term == "Mon, Feb 10 2025") %>%
                       select(N) %>%
                       rename("N of firms" = "N")) %>%
  `class<-`("modelsummary_list")

modelsummary(list("(1) \\textsc{ar}" = mod_ar1,
                  "(2) \\textsc{car}" = mod_car1,
                  "(3) \\textsc{ar}" = mod_ar0,
                  "(4) \\textsc{car}" = mod_car0,
                  "(5) \\textsc{ar}" = mod_ar_diff,
                  "(6) \\textsc{car}" = mod_car_diff),
             # statistic = "[{conf.low}, {conf.high}]",
             notes = paste0("Average \\textsc{ar} and \\textsc{car} to past FCPA targets who did vs did not change name or ownership since FCPA action, ",
                            "per day. Standard errors of the mean reported in parentheses. P-values from a two-tailed test ",
                            "of difference from zero for the average against a standard normal distribution. ",
                            "Estimation window starts 30 days and ends 5 days before FCPA Executive Order. Market models estimated using ",
                            "the LASSO and individual S\\&P 500 constituents as predictors, selected using 15-fold cross validation. ",
                            "Columns 5 and 6 report the difference in means, respectively, between averages in columns 1 and 3, and those in columns 2 and 4.",
                            collapse = ""),
             title = "Heterogeneous effect on past FCPA targets by whether they changed name since FCPA action \\label{tab:het_name_change}",
             stars = c("*" = 0.05), 
             escape = FALSE) %>%
  group_tt(j = list("Changed name" = 2:3,
                    "Did not change name" = 4:5,
                    "Difference-in-means" = 6:7)) %>%
  group_tt(i = list("Pre-event:" = 1,
                    "Post-event:" = 11)) %>%
  style_tt(i = 21, line_color = "white", line_width = 0.1, line = "t") %>%
  style_tt(i = 22, line_color = "black", line_width = 0.05, line = "b") %>%
  theme_tt("resize", width = .9) %>%
  theme_tt("placement", latex_float = "!htbp") %>%
  save_tt("tables/table_M3.html", overwrite = TRUE)

#====# The End #====#